Evidence of an Upper Bound on the Masses of Planets and Its Implications for Giant Planet Formation

I use clustering algorithms from each of three classical clustering paradigms that span the space of modern clustering approaches: connectivity-based hierarchical clustering (Murtagh 1985), centroid-based k-means clustering (Hartigan & Wong 1979), and distribution-based Gaussian-model clustering (Fraley & Raftery 2002; Fraley et al. 2012; Schlaufman 2015). I summarize the key aspects of each algorithm below and provide the details in Appendix B.

The hierarchical clustering algorithm seeks to identify the maximally connected subgraph of a set of objects. It starts with a matrix of dissimilarity between all objects in a set. It then identifies the pair of objects with the smallest dissimilarity, combines them, and updates the dissimilarity of all remaining objects in the set with the minimum dissimilarity to either of the combined objects. It continues this process until only two objects remain. This algorithm requires the specification of a dissimilarity metric, and I use the Euclidean distance between each object.

The k-means clustering algorithm seeks to minimize the total dissimilarity between the ensemble of objects in each cluster and the cluster center. It starts by randomly placing a set of cluster centers in the space defined by the objects to be clustered. It then identifies the closest cluster center to each object and assigns that object to that cluster. The mean value of the ensemble of objects in each cluster is then set to be that cluster's center on the next iteration. The process continues until the movement of cluster centers on each update step is small. This algorithm requires the specification of a dissimilarity metric, and I use the Euclidean distance between each object.

The Gaussian-model clustering algorithm seeks to identify the mixture of Gaussians that has the highest likelihood of having generated the objects to be clustered. It starts by randomly choosing the mean and variance of each Gaussian as well as the contribution of each Gaussian to the mixture. It then calculates the probability that each observation was produced by one component of the mixture. The object is assigned to the component most likely to have generated it and the mean and variance of each component are estimated based on the objects assigned to it. The proportion of each component is the fraction of objects assigned to it. This process is repeated until the change in model parameters is small.

To account for the observational uncertainty in each secondary mass and primary metallicity, I use a Monte Carlo simulation. First, I generate ${10}^{3}+1$ realizations of each secondary mass and primary metallicity from the observational uncertainties on those quantities and then apply all three clustering algorithms to each realization. For each algorithm, I aggregate the classifications across all iterations and assign to each system the most frequently determined classification. The final classification of each system is then the majority classification among the three different algorithms. I plot the result of this calculation in Figure 1. All systems with a secondary mass of ${M}_{2}\lt 4\,{M}_{\mathrm{Jup}}$ are unanimously classified with giant planets, while all systems with a secondary mass of ${M}_{2}\gt 10\,{M}_{\mathrm{Jup}}$ are unanimously classified with low-mass stars. A two-sample Kolmogorov–Smirnov test on the overall sample split at ${M}_{2}=10\,{M}_{\mathrm{Jup}}$ shows that the chance that the two metallicity distributions were produced by the same parent distribution is only $p=4.45\times {10}^{-3}$, or about $3.3\sigma$.¹

Zoom In Zoom Out Reset image size

The clustering analysis described above implicitly makes use of the relative occurrence rates of giant planets, brown dwarfs, and low-mass stars. An alternative approach that does not implicitly use the occurrence rate is to calculate the median metallicity as a function of secondary mass in a kernel of a constant width. Since celestial bodies with $M\gtrsim 10\,{M}_{\mathrm{Jup}}$ formed via core accretion are thought to occur only around the most metal-rich stars (Mordasini et al. 2012), the mass at which the moving median metallicity moves below the lowest moving median metallicity observed in the giant planet region should also correspond to the maximum mass of the objects formed by core accretion.

To account for the observational uncertainty in each secondary mass and primary metallicity, I use a Monte Carlo simulation. First, I generate ${10}^{3}+1$ realizations of each secondary mass and primary metallicity from the observational uncertainties on those quantities and then calculate the moving median metallicity of each sample as a function of mass. I use a kernel of width n = 23 objects, as this is the highest mass resolution that permits the measurement of the median metallicity with a similar precision as for the individual metallicity estimates (${\sigma }_{[\mathrm{Fe}/{\rm{H}}]}\approx 0.05$ dex). I plot this result in Figure 2. Again, the point at which secondaries no longer preferentially orbit metal-rich solar-type dwarf-star primaries occurs at ${M}_{2}\approx 10\,{M}_{\mathrm{Jup}}$. This result is similar for all kernel widths that permit the measurement of the median metallicity with a similar precision as for the individual metallicity estimates.

Zoom In Zoom Out Reset image size

The results in Sections 3.1 and 3.2 show that low-mass secondaries no longer preferentially orbit metal-rich solar-type dwarf stars when ${M}_{2}\gtrsim 10\,{M}_{\mathrm{Jup}}$. Those inferences support the idea that the most massive objects formed around solar-type stars like planets through core accretion have $M\approx 10\,{M}_{\mathrm{Jup}}$. I now use that inference to quantify the range of protoplanetary disk properties that could be expected to reproduce this result using the semi-analytic model of Tanigawa & Tanaka (2016).

Tanigawa & Tanaka (2016) modeled the last stage of giant planet formation (runaway envelope accretion) in a globally evolving protoplanetary disk using the Shakura & Sunyaev (1973) model for disk viscosity and a self-similar solution for the disk's global evolution (e.g., Lynden-Bell & Pringle 1974; Hartmann et al. 1998). The key parameters of the model are as follows.

• 1.  
  The location of planet formation, a_p.
• 2.  
  The disk viscosity, $\nu =\alpha {ch}$, where α is a unitless parameter independent of radius and time, c is the sound speed of disk gas, and h is the disk scale height.
• 3.  
  The initial mass of the disk, ${M}_{\mathrm{disk}}$, is assumed to scale linearly with a host star mass, M_*, and overall normalization, ${f}_{{\rm{\Sigma }}}$. Here ${f}_{{\rm{\Sigma }}}=1$ corresponds to the minimum-mass solar nebula of Hayashi et al. (1985).
• 4.  
  The characteristic exponential disk depletion time, ${\tau }_{\mathrm{dep}}$.

The model is informed by the latest numerical results showing that gaps cleared by giant planets are much shallower than previously thought (e.g., Duffell & MacFadyen 2013; Fung et al. 2014). The model does not account for possible photoevaporation of disk gas due to far-ultraviolet radiation (e.g., Gorti et al. 2009), as the rate of mass loss expected is still uncertain by more than an order of magnitude (e.g., Alexander et al. 2014).

I have implemented the Tanigawa & Tanaka (2016) model and used it to conduct a parameter study of the relationship between the maximum attainable mass of a giant planet and the model parameters a_p, α, M_*, ${f}_{{\rm{\Sigma }}}$, and ${\tau }_{\mathrm{dep}}$. If not otherwise indicated, I assume that a_p = 1 au, ${M}_{* }=1\,{M}_{\odot }$, ${\tau }_{\mathrm{dep}}\,=3\times {10}^{6}\,\mathrm{year}$, the initial disk radius is ${R}_{0}=200\,\mathrm{au}$, and the seed embryo mass is ${M}_{\mathrm{embryo}}=10\,{M}_{\oplus }$. I plot the result in Figure 3. According to the model, a maximum giant planet mass in the range of $4\,{M}_{\mathrm{Jup}}\lesssim M\lesssim 10\,{M}_{\mathrm{Jup}}$ can be accommodated only in a narrow corridor of the α–${f}_{{\rm{\Sigma }}}$ plane from $\alpha \sim {10}^{-4}$ and ${f}_{{\rm{\Sigma }}}\sim {10}^{0.5}$ to $\alpha \sim {10}^{-2}$ and ${f}_{{\rm{\Sigma }}}\sim {10}^{-1}$. This result is only weakly dependent on the mass of the primary star, the location of planet formation, or the disk depletion time. These results are unique in that they indirectly constrain the viscosities and masses of long since dissipated disks during their epoch of giant planet formation.

Zoom In Zoom Out Reset image size

The results presented in Section 3 indicate that a $M\approx 10\,{M}_{\mathrm{Jup}}$ maximum mass for objects formed by core accretion can only be accommodated in a narrow band in the Shakura & Sunyaev (1973) viscosity parameter α–minimum-mass solar nebular scaling factor ${f}_{{\rm{\Sigma }}}$ plane from $\alpha \sim {10}^{-4}$ and ${f}_{{\rm{\Sigma }}}\sim {10}^{0.5}$ to $\alpha \sim {10}^{-2}$ and ${f}_{{\rm{\Sigma }}}\sim {10}^{-1}$. Observations of T Tauri stars have been used to suggest that protoplanetary disks can be reproduced by models with $\alpha \sim {10}^{-2}$ (e.g., Hartmann et al. 1998). At the same time, classical models of giant planet formation though core accretion almost always invoke disks several times more massive than that of the minimum-mass solar nebula—or ${f}_{{\rm{\Sigma }}}\gtrsim 3$—to ensure core formation before disk dissipation (e.g., Pollack et al. 1996; Ikoma et al. 2000; Alibert et al. 2005; Hubickyj et al. 2005; Lissauer et al. 2009; Movshovitz et al. 2010; D'Angelo et al. 2014). In a disk with $\alpha ={10}^{-2}$ and ${f}_{{\rm{\Sigma }}}=3$, the Tanigawa & Tanaka (2016) model predicts a giant planet maximum mass of $M\approx 60\,{M}_{\mathrm{Jup}}$. According to these results, either but not both, of these traditional assumptions could be made during the runaway gas accretion stage of giant planet formation.

If the $M\approx 10\,{M}_{\mathrm{Jup}}$ threshold is the result of low-viscosity or low-mass disks, then it would provide a solution to the long-standing Type I and Type II migration problems. Type I migration occurs when planetary embryos with ${M}_{\mathrm{embryo}}\,\sim 1\,{M}_{\oplus }$ that are not massive enough to open a gap in their parent protoplanetary disks lose angular momentum and move inward due to net negative torques from the disk. Though the processes responsible are complex and no model is yet completely satisfactory, the simplest models suggest that the timescale for migration is very fast (e.g., Ward 1997; Tanaka et al. 2002). Exoplanet population synthesis calculations assuming fast Type I migration are unable to reproduce the existence of planets with $M\gtrsim 10\,{M}_{\oplus }$ and $P\gtrsim 10$ days, so Type I migration must be slower than the simplest calculations initially suggested (e.g., Ida & Lin 2008; Schlaufman et al. 2009). While many solutions have been proposed,² it has been shown that Type I migration in disks with $\alpha \sim {10}^{-4}$ is sufficiently slow that planetary embryos will survive and potentially grow into giant planets (e.g., Rafikov 2002; Li et al. 2009). Similarly, in disks with $\alpha \sim {10}^{-2}$, the Type I migration rate is thought to be linearly dependent on disk surface density (e.g., Baruteau et al. 2016). A minimum-mass solar nebular scaling factor of ${f}_{{\rm{\Sigma }}}\sim {10}^{-1}$ would reduce the Type I migration rate by an order of magnitude, exactly the amount needed to reproduce the observed distribution of exoplanet properties in population synthesis calculations (e.g., Ida & Lin 2008; Schlaufman et al. 2009).

Low-viscosity or low-mass disks would also solve the Type II migration problem. Type II migration occurs when planets that are massive enough to open a gap in their parent protoplanetary disks are carried into the close proximity of their host star on the disk viscous time (Lin et al. 1996), possibly precluding the existence of long-period giant planets. In disks where the mass of the disk locally dominates the mass of the planet, the Type II migration rate is inversely proportional to the viscosity (e.g., Lin et al. 1996; Tanigawa & Tanaka 2016). Similarly, in a low-mass disk where the mass of the planet locally dominates the mass of the disk, the Type II migration rate is inversely proportional to the disk mass (e.g., Duffell et al. 2014). Consequently, low-viscosity or low-mass disks promote the existence of long-period giant planets.

Though these results were derived strictly for solar-type dwarf stars with short-period transiting companions, the analysis in Section 3.3 suggests that the inference about planet formation is not strongly dependent on host-star mass, planet-formation location, or disk dissipation timescale. Indeed, similar inferences have been made using the much-more distant solar system giant planets Jupiter and Saturn by Lissauer et al. (2009) and Tanigawa & Tanaka (2016). It is also reassuring that this upper limit of $10\,{M}_{\mathrm{Jup}}$ on the mass of planets within 1 au of their host stars is consistent with the most advanced exoplanet population synthesis calculations currently available (Ida et al. 2013). Finally, support for low-viscosity or low-mass disks as the origin of planetary systems has been found in many places from meteoritics in the solar system to super-Earths and mini-Neptunes in the Kepler field (e.g., Desch et al. 2017; Lee et al. 2017).

This study focused on giant planets, brown dwarfs, and low-mass stars with Doppler-inferred masses transiting solar-type dwarf stars with homogeneously estimated stellar parameters. It showed that low-mass secondaries with ${M}_{2}\lt 10\,{M}_{\mathrm{Jup}}$ preferentially orbit metal-rich stars, while secondaries with ${M}_{2}\gt 10\,{M}_{\mathrm{Jup}}$ occur independent of host star metallicity. The difference between the metallicity distributions of objects in my sample with $M\gt 10\,{M}_{\mathrm{Jup}}$ and $M\lt 10\,{M}_{\mathrm{Jup}}$ is significant at more than $3.3\sigma$, and this result cannot be attributed solely to the possible change in the occurrence rate of such bodies at $M\approx 4\,{M}_{\mathrm{Jup}}$. Since a preference for metal-rich host stars is thought to indicate formation through core accretion, these results suggest that the most massive objects that can be formed like planets through core accretion have $M\approx 10\,{M}_{\mathrm{Jup}}$.

Unlike all previous studies, these results are unaffected by the Doppler $\sin i$ degeneracy, the presence of possible false positives in the analysis sample, or the potentially weak or nonexistent giant planet–host star metallicity correlation for giant stars (e.g., Sadakane et al. 2005; Schuler et al. 2005; Hekker & Meléndez 2007; Pasquini et al. 2007; Takeda et al. 2008; Ghezzi et al. 2010; Maldonado et al. 2013; Mortier et al. 2013a; Jofré et al. 2015; Reffert et al. 2015). Because of its large sample and use of M instead of $M\sin i$, it has the best resolution in mass of all studies to date. All of these properties have enabled the identification of a difference between the metallicity distributions of objects with $M\gt 10\,{M}_{\mathrm{Jup}}$ and $M\lt 10\,{M}_{\mathrm{Jup}}$ at a higher statistical significance than any previous study of solar-type dwarfs. It is also the first study to couple its result to a detailed model of the runaway gas accretion stage of giant planet formation to constrain the properties of protoplanetary disks that are known to have formed giant planets.

The first studies seeking to compare giant planets and brown dwarfs to explore their origin all focused on the orbital properties of brown dwarfs. Grether & Lineweaver (2006) examined the mass distribution of companions to solar-type stars with $P\lt 5\,\mathrm{year}$ and determined that the minimum of the distribution occurs at $M\sin i={31}_{-18}^{+25}\,{M}_{\mathrm{Jup}}$. Because there was no such minimum in the mass spectrum of free-floating brown dwarfs, they speculated that migration and not a characteristic mass scale of fragmentation in protoplanetary disks was responsible for the apparent minimum.

Sahlmann et al. (2011) used the combination of the Doppler and astrometric techniques to remove the $\sin i$ degeneracy to confirm the previous result, finding that low-mass brown dwarfs with $M\lesssim 25\,{M}_{\mathrm{Jup}}$ cleanly separated from high-mass brown dwarfs with $M\gtrsim 45\,{M}_{\mathrm{Jup}}$. In contrast to Grether & Lineweaver (2006), they argued that the low-mass brown dwarfs may be massive planets formed by core accretion. They also showed that samples of brown dwarf candidates discovered solely by the Doppler technique with $M\sin i\,\gtrsim 45\,{M}_{\mathrm{Jup}}$ are significantly contaminated by inclined low-mass stars.³

Ma & Ge (2014) extended these results by demonstrating that brown dwarfs with $M\sin i\lesssim 42.5\,{M}_{\mathrm{Jup}}$ and $M\sin i\,\gtrsim 42.5\,{M}_{\mathrm{Jup}}$ are unlikely to have the same eccentricity distribution at the 2.1σ level. They suggested that the low-mass brown dwarfs formed by disk fragmentation, while the high-mass brown dwarfs formed by fragmentation at the core scale. They also showed that the giant planet and brown dwarf host stars were unlikely to have the same metallicity distribution at the 3.3σ level. Selection effects in the giant planet–host sample and inclined low-mass star contamination in the brown dwarf sample likely affect these conclusions to some degree though.

More recent efforts to investigate the differences between massive planets and brown dwarfs have focused on host star metallicity and chemical abundance. Using a sample of solar-type dwarf stars unaffected by selection effects or inclined low-mass star contamination, Mata Sánchez et al. (2014) found that the metallicities and α-element abundances of brown dwarf hosts are lower than the metallicities of giant planet–host stars at the 1.7σ and 1.8σ levels.

Maldonado & Villaver (2017) presented a similar analysis to Ma & Ge (2014) that used homogeneously derived stellar parameters to show that low-mass brown dwarf hosts are more metal-rich than high-mass brown dwarf hosts at the 2.4σ level. However, the difference was not significant if only dwarf host stars were included or if a random distribution of inclination was assumed. The fact that the significance of their result is reduced when accounting for a random distribution of inclination suggests that low-mass stars on low-inclination orbits contaminated their sample. In addition, their compete sample of brown dwarf hosts was significantly more metal-poor than their sample of giant planet–host stars. At the same time, their low-mass brown dwarfs were consistent with having the same metallicity distribution as the giant planet–host stars, though their α-element distributions differed at the 2.4σ level.

Santos et al. (2017) studied the minimum-mass and metallicity distributions of mostly Doppler-discovered objects with $M\sin i\lt 15\,{M}_{\mathrm{Jup}}$, $10\,\mathrm{days}\lt P\lt 1825\,\mathrm{days}$, and with homogeneous stellar parameters available in the SWEET-Cat database. Using both multivariate Gaussian modeling and Kolmogorov–Smirnov tests, they argued that planets below and above $M\sin i=4\,{M}_{\mathrm{Jup}}$ have distinct metallicity distributions and therefore distinct formation channels. In particular, they suggested that the giant planets with $M\sin i\lesssim 4\,{M}_{\mathrm{Jup}}$ formed via core accretion while the giant planets with $M\sin i\gtrsim 4\,{M}_{\mathrm{Jup}}$ formed via gravitational instability. Their multivariate Gaussian inference is mostly due to the apparent cutoff in the giant planet occurrence distribution at $M\sin i\approx 4\,{M}_{\mathrm{Jup}}$. In their Kolmogorov–Smirnov tests, the distributions were distinct at the 3.3σ significance substantially because of the inclusion of giant stars with $M\gtrsim 1.5\,{M}_{\odot }$. The fact that restricting the sample to solar-type dwarf stars reduced the significance of their inference to 1.6σ indicates that either false positives in the giant star sample or the possibly diminished giant planet occurrence–host star metallicity effect for giant stars may play a role in the significant difference they identified in their sample including both dwarf and giant host stars. A critical mass of $M\sin i=4\,{M}_{\mathrm{Jup}}$ is also inconsistent with the heavy element enrichment observed in HAT-P-20 b and HAT-P-2 b. Nevertheless, my study confirms the hint identified in Santos et al. (2017) and presents the strongest evidence to date that the maximum mass of objects formed via core accretion is in the range of $4\,{M}_{\mathrm{Jup}}\lesssim M\lesssim 10\,{M}_{\mathrm{Jup}}$.

This analysis has shown that objects with $M\lesssim 10\,{M}_{\mathrm{Jup}}$ preferentially orbit metal-rich solar-type dwarf stars. That property has been suggested to be a natural outcome of the core accretion model of giant planet formation. On the other hand, objects with $M\gtrsim 10\,{M}_{\mathrm{Jup}}$ orbit stars spanning the whole range of the thin disk metallicity distribution. That property is shared with low-mass stars, which must have formed by some sort of gravitational instability.

I propose that planets be defined as objects that orbit stars or stellar remnants and have a true mass below the threshold at which low-mass companions no longer preferentially orbit metal-rich solar-type dwarf stars. The analyses presented in Sections 3.1 and 3.2 suggest this threshold is at $M\approx 10\,{M}_{\mathrm{Jup}}$, but future data may revise this estimate. This definition has the same form as the working definition adopted by the IAU Working Group on Extrasolar Planets, but the maximum mass is now referenced to the maximum mass of objects formed by core accretion. This definition does not require the specification of an arbitrarily amount of deuterium burning, and it does not depend on an object's metallicity. Additionally, it is independent of the uncertain internal structure of objects at this mass scale. Furthermore, I propose that substellar objects with true masses above the threshold be defined as brown dwarfs, regardless of where they are located.

The results presented in Sections 3.1 and 3.2 lead to the prediction that planets with $1\,{M}_{\mathrm{Jup}}\lesssim M\lesssim 10\,{M}_{\mathrm{Jup}}$ should preferentially be found around metal-rich solar-type dwarf stars. In contrast, brown dwarfs with $10\,{M}_{\mathrm{Jup}}\lesssim M\lesssim 80\,{M}_{\mathrm{Jup}}$ should be found around solar-type dwarf stars that span the metallicity range of the Milky Way's thin disk. This prediction will soon be tested at an unprecedented scale. By the end of its planned five-year mission, the European Space Agency's Gaia satellite is expected to discover astrometrically and characterize 21,000 ± 6000 giant planets and brown dwarfs with $1\,{M}_{\mathrm{Jup}}\lt M\lt 15\,{M}_{\mathrm{Jup}}$ and $1\,\mathrm{au}\lesssim a\lesssim 5\,\mathrm{au}$ (Perryman et al. 2014). This sample of giant planet and brown dwarf mass measurements will be nearly 100 times larger than the sample currently available and extend to the more common intermediate-period objects lacking in my sample. For these reasons, Gaia will provide the definitive test of this prediction that planets with $1\,{M}_{\mathrm{Jup}}\lesssim M\lesssim 10\,{M}_{\mathrm{Jup}}$ should preferentially be found around metal-rich solar-type dwarf stars, while brown dwarfs with $10\,{M}_{\mathrm{Jup}}\lesssim M\lesssim 80\,{M}_{\mathrm{Jup}}$ should be found around solar-type dwarf stars that span the metallicity range of the Milky Way's thin disk.

I use the hclust function in R for hierarchical clustering (R Core Team 2017). That implementation follows the algorithm of Murtagh (1985).

• 1.  
  Calculate the set of $n(n-1)/2$ dissimilarities between n objects.
• 2.  
  Identify the smallest dissimilarity d_ik.
• 3.  
  Combine objects i and k. That is, replace the two objects with a new object, $i\bigcup k$, and update all other dissimilarities such that for all objects $j\ne i,k$ the dissimilarity ${d}_{i\bigcup k,j}=\min ({d}_{i,j},{d}_{k,j})$. Delete the dissimilarities ${d}_{i,j}$ and ${d}_{k,j}$.
• 4.  
  Return to Step 2 while the number of objects remaining is greater than two.

This algorithm requires the specification of a dissimilarity metric, and I use the Euclidean distance between each object.

I use the kmeans function in R for k-means clustering (R Core Team 2017). That implementation follows the algorithm of Hartigan & Wong (1979).

• 1.  
  Pick k objects from the n objects to be clustered.
• 2.  
  Assign each object to a cluster by placing each it in the cluster, k_i, that has the centroid with the smallest dissimilarity.
• 3.  
  Recalculate the centroid of each cluster.
• 4.  
  Return to Step 2, unless changes in the cluster centroids are small.

This algorithm requires the specification of a dissimilarity metric, and I use the Euclidean distance between each object.

I use the mclust package in R for Gaussian-model clustering (Fraley et al. 2012; R Core Team 2017). The algorithm is considerably more complex than either the hierarchical or k-means algorithms, so the interested reader should see Fraley & Raftery (2002) for the details of the algorithm.